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Abstract 

We present the exact and precise (~0.1%) numerical solution of the QCD evolu- 
tion equations for the parton distributions in a wide range of Q and x using Monte 
Carlo (MC) method, which relies on the so-called Markovian algorithm. We point 
out certain advantages of such a method with respect to the existing non-MC meth- 
ods. We also formulate a challenge of constructing non-Mar kovian MC algorithm 
for the evolution equations for the initial state QCD radiation with tagging the 
type and x of the exiting parton. This seems to be within the reach of the presently 
available computer CPUs and the sophistication of the MC techniques. 
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1 Introduction 



It is commonly known that the so called evolution equations of the quark and gluon 
distributions in the hadron, derived in QED and QCD using the renormalization group or 
diagrammatic techniques pj, can be interpreted probabilistically as a Markovian process. 
This process consists of the random steps forward in the logarithm of the energy scale, 
t = InQ, and the coUinear decay of partons. Such a Markovian process can be readily used 
as a basis of the Monte Carlo (MC) algorithm producing events in the so called Parton 
Shower MC (PSMC), see for example ref. |2j. The Markovian MC provides, in principle, 
an exact solution of the evolution equations. This possibility was not exploited in practice 
in the past, mainly because alternative non-MC numerical methods and programs solving 
the QCD evolution equations on the finite grid of points in the space of t and parton 
energy x are much faster than the MC method. Typical non-MC program of this type is 
QCDnumie of ref. [3]. For comparisons with other codes see ref. P|. 

In this work we test practical feasibility of applying the MC Markovian method to solve 
the QCD evolution equations exactly and precisely (per- mill level), exploiting computer 
CPU power available today, in a wide range of InQ and x. 

Apart from pure exploratory aspect, we believe that MC method has certain advan- 
tages over non-MC methods, see discussion below. 

It is also commonly known that the basic Markovian algorithm cannot be used in the 
PSMC to model the initial state radiation (ISR), because of its extremely poor efficiency. 
In the latter part of this work we formulate a challenge of constructing non-Markovian type 
of MC algorithm for ISR PSMC, in which the exact evolution of the parton distributions 
is a built-in feature, like in the Markovian MC case^. We claim that with present CPU 
power and sophistication of the MC techniques such a scenario may be feasible. With the 
advent of the future high quality and high statistics experimental data coming from LHC 
and HERA, such a solution is worth to consider. 

The layout of the paper is the following: In Section 2 we summarize basic concepts on 
the Markovian process and its application to solution of the QCD evolution equations. In 
Section 3 we present exact high precision Markovian MC solution of the QCD evolution 
equations. Section 4 discusses the results and perspectives. 

2 The formalism 

Let us briefly review basic definitions and concepts concerning the evolution equations 
and the Markovian process. The generic evolution equations read 



K 

where t = \nQ and / and K are discrete indexes numbering states of a certain physical 
system. The extension to continuous sets of states is straightforward, see below. The 

^This is contrary to most of the existing MC PSMCs, in which evolved parton distributions come from 
a table obtained using an independent non-MC program. 
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necessary and sufficient condition for this equation to have its representation as a proba- 
bihstic Markovian process in the evolution time t, with a properly normalized probability 
of a forward step, is the following: 



where Ri{t) > is the decay rate of the I-th parton. The evolution equation 

dtDj{t) + Ri{t)Di{t) = J2 'yiK{t)DK{t), 

K^I 

can be easily brought to a homogeneous form 

K^I 

^i{t,to) = I dti Riih), 

J to 

which can be turned into an integral equation 

Dj{t) = e-*^(*'*°)L',(io) + fdh e-*^(*'*^) ^^^(^i) ^i^(^i) 



•J I 



KJ 



Tkj, for K^J, 
0, for K = J, 



and finally can be solved by means of multiple iteration: 

oo r 
i^xW = e-*-(*'*°)L>^(io)+ J] En/ dUe{U-t,_^ 

n=l Ko-Kn-i i=l '-'^*o 
n 



i=l 



(2) 



(3) 



(4) 



(5) 



(6) 



-Dxo(^o), 



where K = is understood, for the brevity of the notation. The above series of integrals 
with positively defined integrands can be interpreted in terms of a random Markovian 
process starting at to and continuing until t with the following, properly normalized, 
transition probability of a single Markovian random step forward, (ij-i, i^^i-i) {U, Ki), 



/•oo 

/ dtj y^^Lj{ti,Ki\ti_i,Ki_i) = 1. 
Jti-i 



(7) 
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Each term in the sum in eq. © can be expressed as a product of the single step proba- 
bihties of eq. ((Zj) as follows: 

DK{t) =e-*-(*'*°)D^(to) 

n=l Ko...K„.ij=l-^^0 j=i 

The Markovian process does not run forever - it is stopped by the so called stopping rule, 
which in our case is tn+i > t. The variable tn+i is added in eq. (jH)) by means of inserting 
an extra integration with the help of an identity 



oo 



J2uj{tn+l,Kn+l\tn,K^). (9) 

^ Kn + l 

The resulting solution 

DKit)= I dti^uj{ti,Ki\to,Ko)DK{t^) 
tit ^1 

OO « Up 

+ ^ ^ / dtn+1 ^ Uj{tn+l,Kn+l\tn,Kn)Y\ / ^tj (10) 
n=l Ko...K„-i Kn+i ^=\<t 

n 

1=1 



represents clearly a Markovian process, in which, starting from a point (to, -^o) distributed 
according to DKoito), we generate step by step the points (tj, Ki),i = 1, 2, . . . , n + 1, using 
transition probability u, until the stopping rule tn+i > t acts. In such a case the point 
[tn+i, Kn+i) is trashed and the previous point (tn,Kn) is kept as the last one. Eq. (fTIH) 
tells us that the points (t„, K^) obtained in the above Markovian process are distributed 
exactly according to the solution of the original evolution equations^ of eq. ((T)). 

Two technical point have to be clarified before we apply the above solution to the 
QCD parton distribution evolution equations: one concerns the t dependence of the 5" 
matrix and another one concerns extension of the above formalism to continuous indexes 
in the T matrix. 

The transition matrix T, in the case of the QCD parton distributions, includes the 
running coupling constant a(t) 

nt)KL = ait) Pkl, (11) 



^The above MC solution is exact and unbiased contrary to other non-MC methods which have inherent 
biases due to choice of finite grid decomposition into polynomials and other technical artifacts present 
there. The main disadvantage of the MC method is its slowness. 
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which depends rather strongly on t at low t = \nQ. In the above we have explicitly 
assumed that the kernel Prl is independent of t, which is true in the leading-logarithmic 
(LL) case and is not true in the next-to-leading-logarithmic (NLL) case. However, even 
in the NLL case this is a very good approximation, which can be easily corrected in the 
MC calculation using an extra weight (and event rejection). So, we can safely assume 
the validity of eq. |TT| without any loss of generality, even beyond the LL level. If the 
above factorization (fTT|) is true, then we may employ the standard trick which eliminates 
t-dependence from T completely. This is realized by introducing a new evolution time 
variable 

1 /■* , r ^ dt a(to) 
a{to) Jt^ dT a[t) 

With this choice we have 

drDi{T) = J2 'PiKDKir), TiK = ^'^iKit) = a{to)PKL, (13) 



K 



where T does not depend on the new evolution time variable r anymore, because a{t) in 
y is effectively replaced by a(to)- The rest of the formalism following eq. is the same, 
provided we substitute t — > r and 7 T everywhere. 

Concerning the continuous indexes, indeed in the QCD case we deal with the composite 
index K = {k,x}, where k is the parton type; let it be /c = G,q,q in the simple case 
of gluon and one type of quark, and x, which is a fraction of the energy of the primary 
(proton) beam particle carried by the parton - the usual definition. The standard parton 
distribution evolution equations read 

1 

drXiDk^{T,Xi) [ dXo "^^'^"^ XoDko{T,Xo) (14) 

from which we deduce the following substitution rule 

while for the matrix elements, excluding IR region Xi = Xi^i from the consideration, we 
have the explicit expression 

7k,k,_, = —Pk.k^.^i B 1 , 16 

TT Xf_^ \Xi-iJ \ Xi-iJ 

where = {ki,Xi}. 

Another two points require now clarification: the meaning of the diagonal part K = J 
in the the off- diagonal part and the validity of the condition of eq. necessary 
for the feasibility of the Markovianization. 
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The problem of defining off-diagonal elements in T'^ j coincides with the problem of 
defining an infrared (IR) cut-off and is dealt with in a standard way by introducing a 
dimensionless small parameter e — 0, which defines/separates real gluon emission by 
requiring < 1 — e. The off-diagonal (real emission) part is defined in the usual 

way^ 

= e(i-.-|)^|p,^.^(|) +o(.)(i-4..,). (17) 

The condition of eq. (j2I), J2k ^KjKi = 0, coincides in fact with the energy conservation 
sum rule, which reads as follows 

Kj kj 

where Zji = Xj/xi, and the nice feature is that depends only on the parton type 
ki and not on Xi (unless we introduce more complicated IR cut). At this point it is 
understandable why eq. ()14|) was written for the energy distributions xD{x) and not for 
the parton distributions Dk{x) themselves. In the latter case the condition of eq. (j21) is 
not generally valid^ and the Markovian MC calculation is not possible. The significance 
of the energy sum rules in the evolution equations was underlined in many works, see for 
instance refs. IH] . 

Collecting all the above material we can write the normalized transition probability 
for the actual MC Markovian calculation as follows 

uj{Ti, ki, x,|ri_i, Xi_i) = e(r, - r,_i) g-^^^^^-i)-^'-! (19) 

where the necessary ingredients are provided by eqs. (fT7j) and (fTHj) . The LL QCD evolution 
kernels can be found in any QCD textbook and we skip details of the actual formulas 
entering the above equation. 

3 Numerical results 

In the following we present numerical results obtained for the distributions of three light 
quarks and gluons evolved from Q = IGeV up to ITeV in the range x G (10~^, 1). 
The starting parton distribution at Qo = IGeV is the proton distribution defined in a 

•^The 0{e) term is related to the fact that for the sake of simphcity we apply the Q{xj/xi — 1 + e) 
factor also to the ki kj part, which is not IR singular. 
'*It is, however, valid for the non-singlet component. 
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Figure 1: The upper plot shows gluon distribution xDcix.Qi) evolved from Qq = iGeV to 
Qi = 10, 100, lOOOGeV obtained from QCDnuml6 and EvolMCl, while the lower plot shows 



their ratio. The horizontal axis is log 



10 



See the text for details. 



conventional way: 



xDg 


'x) 


= 1.9083594473 


.a;-o.2(i_3,)5.o 


XDg 


'x) 


= 0.5 ■ xDsc!,{x) 


+ xD2uix), 


XDg 


[x) 


= 0.5 • xDsei,{x) 


+ xDd{x), 


xDsea, 


'x) 


= 0.6733449216 


.x-^-\i-xy-' 


xD2u 


[x) 


= 2.1875000000 


■x'-\l-xf'>, 


xDd 


[x) 


= 1.2304687500 


■x'-\l-x)^-°. 



(20) 



We use in this exercise the LL version of the evolution kernels |TJ [7j . Extension to the 
NLL level requires to use the NLL kernels |HI. The rest of the MC algorithm/program is 
the same as in the LL case. Such an extension will be presented elsewhere. 

In figs. we show numerical results for the gluon and quark distributions evolved 
from Q = IGeV to higher scales Q = 10, 100, lOOOGeV, obtained from our Markovian 
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Figure 2: The upper plot shows the singlet quark distribution D 



X 



evolved from Qq 



iGeV to Qi = 10, 100, lOOOGeV obtained from QCDnumie and EvolMCl, while the lower plot 



shows their ratio. The horizontal axis is log 



10 



x] 



See the text for details. 



Monte Carlo EvolMCl^. In the plots we also show the results from the standard program 
QCDnuml6 evolving parton distributions using a finite grid of points in the x and Q space. 
As we see, these two calculations agree within 0.2% for gluons and even better for the 
quark singlet. This is clearly seen from the plotted ratio of the MC result and the 
result of QCDnuml6. The biggest discrepancy is in the region close to x = 1, where all 
parton distributions are extremely small. The origin of 0.2% discrepancy for gluons is 
not identified unambiguously. Variation of the grid parameters change the results of 
QCDnumie quite substantially, while the variation of the IR cut e and of the minimum 
X parameters in the MC run has no visible influence on the MC results. We tend to 
conclude that the difference between the MC and QCDnuml6 is due to a numerical bias in 
the QCDnuml6 program, although more work is required to reach a firm conclusion. In the 
MC we generated 1.2-10^ events (96nh CERN units). It would be possible to increase MC 
statistics by a factor of 10, obtaining a sub per-mill statistical errors in the MC results, 

^The actual MC program is rather compact, it consists of about 3k lines of the C++ code. 
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if necessary. 



4 Discussion 




logio(x) logio(x) 

Figure 3: Average multiplicities of all emissions (left plot) and of the transitions G — > g and 
g — > G only (right plot) as a function of log^Q(x). Results are for the proton "showering" up 
to ITeV scale according to Markovian MC with the IR cut e = 10~^. 

The main lesson from the above numerical exercises is that with the modern computing 
power the precision of 10~^ in the solutions of the QCD evolution equations is within 
the reach of the MC method for most of the a;-range. MC will always be slower than 
the non-MC methods, however, it has no biases related to the finite grid, the use of 
quadratures, decomposition into finite series of polynomials, accumulation of the rounding 
errors, etc. MC method can therefore be used to test any other non-MC tool for the QCD 
evolution. Let us stress that in the MC method there is also no need to split PDFs into 
a singlet component and several types of non-singlet components, each of them evolved 
separately and combined in the very end of the calculation. This simplifies the calculation, 
and the MC Markovian method is clearly well suited for the evolution of PDFs with 
many components, for example QCD+QED, SUSY theory or QCD+electroweak theory 
at very large energies, modeling showers of extremely high energetic cosmic rays in the 
intergalactic space (see ref. 0), etc. 

Let us now elaborate more on the issues related to the construction of the efficient 
PSMC for the QCD initial state radiation. As it is well known, the cross section of 
the hard scattering, especially in the presence of narrow resonances, discriminates very 
strongly on the parton type and x of the parton exiting PSMC and entering the hard 
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scattering. The Markovian PSMC does not provide for any control (tagging) of the flavor 
and X of the exiting parton, consequently, the overall efficiency of such a MC for the 
initial state would be drastically bad. The solution adopted in the most of the existing 
ISR QCD PSMCs is the so-called backward parton shower of ref. JU], in which one starts 
with the parton next to the hard process, defining its type and energy x, and then the 
Markovian random walk continues toward the lower Q and higher x, until the hadron 
mass scale is reached. This beautiful algorithm requires, however, the a priori knowledge 
of the parton distribution functions (PDF) at all intermediate Q scales. These PDFs 
are, therefore, tabulated using results of the non-MC evolution programs, and are also 
fitted to the deep inelastic scattering (DIS) data. In other words, the evolution of PDFs 
is not a built-in feature of the standard ISR PSMC. The two most important reasons 
for that are the following: (i) one is purely technical - the low efficiency of the standard 
Markovian ISR PSMC and another one is (ii) the convenience in porting PDFs from DIS 
data to hadron-hadron scattering. We notice that both of these reasons seem nowadays 
to fade away. For (i), the CPU power of the present computer systems is bigger by the 
factor ~ 10^ than 15 years ago, when the problem of the efficiency of the PSMC for ISR 
was considered for the first time. That widens considerably the range of the practically 
realizable MC algorithms. In addition, there exist now working examples of the efficient 
non-Markovian algorithms with built-in evolution (albeit in the soft approximation) for 
the pure multi-bremsstrahlung process in the abelian case (QED), see refs. |1H IT^. which 
can be treated as a partial solution of the problem of tagging x. For (ii), PDFs in the 
future LHC experiments will be constrained not only by DIS but also by the LHC data {W 
and Z production) themselves. One should, therefore, consider the genuine MC solution 
for the QCD evolution as an integral part of the PSMC, with the aim of fitting the 
LHC and DIS data simultaneously. The obvious advantage of such a scenario is that the 
detector effects can be removed from the data more efficiently. No doubt, to realize it in 
practice will be technically rather challenging. Nevertheless, we claim that it is feasible 
and at present it makes sense to consider it seriously. 

In view of the above discussion let us formulate the challenge in the construction of the 
MC event generators: It would be desirable to construct a MC program/algorithm, which 
evolves the parton distribution functions using the QCD NLL MS kernels from the PDFs 
at low Qmin to higher Qmax, in which we are able to fix beforehand the type (fiavour) and 
the energy fraction x of the parton at the scale Qmax- The a priori knowledge of the PDFs 
at all intermediate Qs should not be required/used; the input from the perturbative QCD 
would be the expressions for the kernels. The PDFs at all intermediate Qs would be the 
result from the MC. 

Additional remarks: (i) We do not assume that the solution is based on the Markovian 
process. Any, even brute force solution, but feasible with the present CPU power, is 
worth consideration, (ii) The above challenge concerns primarily the {Q,x) space. The 
construction of a full scale parton shower MC on top of the above two-dimensional MC 
we treat as a separate (important) issue, (iii) In our opinion, solving the above technical 
challenge would open new interesting avenues in the realm of the MC event generators 
for DIS and hadron-hadron scattering. 
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For the moment we do not offer any solution of our own. We can see one hint, which 
comes from the exercises we have done using the presented MC of the Markovian type. 
In the left hand side plot in Fig. |21 we show the average number of all "emissions" in 
the Markovian process on the way from Q =lGeV to Q =lTeV, for the proton case, 
as a function of the final x at Q =lTeV. The IR cut-off e = 10~^ is used. Due to a 
sizable value of 05 and large phase space the average number of emissions is about 20 
(for a single proton beam). However, the average number of the transitions q,q G and 
G ^ q,q, shown in the right hand side plot in the same figure, is much smaller - it is in 
the range from 0.5 to 1. This feature suggests the brute force solutions in which q,q G 
and G ^ q,q transition are pretabulated (solving the problem of tagging the parton 
type), while the pure bremsstrahlung transition G ^ G and q ^ q are treated using the 
algorithm of the QED (with the constrained final x). However, we hope that more elegant 
and efficient solution can be found. Finally, we would like to refer the reader to extensive 
literature on the subject of the QCD evolution and parton shower MC algorithms, which 
can be found for instance in refs. [71 113j. 
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